
\chapter{矩阵本征值}

矩阵本征方程为
\begin{equation}
A \vec{x} = \lambda \vec{x},
\end{equation}
我们只讨论$A$为实对称阵的情况。
哈密顿量是厄米的，实数的哈密顿量矩阵是实对称阵，它的本征值就是本征能量，本征函数就是体系本征态。

\section{Jacobi方法}
数学家们已经证明了，任意实对称阵$A$都存在正交变换$S$，使得$S^\top A S$为对角阵。
所以对角化任务变成了：构造正交变换$S$，使得$S^\top A S$是对角阵。

Jacobi方法构造一系列旋转，最终使得$A$近似地对角化，完成目标。
在p-q平面上逆时针旋转$\phi$角对应的变换矩阵为
\begin{equation}
P_{pq} = \begin{bmatrix}
 1 \\
  & \ddots \\
  &     & c     &     & -s \\
  &     &       & \ddots \\
  &     & s    &         & c \\
  &     &       &        &       & \ddots \\
  &     &       &        &       &       & 1
\end{bmatrix}
\end{equation}
其中$c=\cos\phi,-s,s=\sin\phi,c$分别处在第$p,q$行，第$p,q$列，矩阵的其余部分与单位矩阵相同。
注意$P_{pq}$是一个矩阵，不是带下标的矩阵元。
在旋转变换以后，矩阵A变为
\begin{equation}
A' = P^\top_{pq} A P_{pq}
= \begin{bmatrix}
        &        &  a'_{0p} &       & a'_{0q} &         &   \\
        &        &   \vdots &       &  \vdots &         &   \\
a'_{p0} & \cdots & a'_{pp}  & \cdots& a'_{pq} &\cdots   & a'_{p,n-1} \\
        &       & \vdots    &       & \vdots  &     \\
a'_{q0} & \cdots & a'_{qp}  & \cdots& a'_{qq} &\cdots   & a'_{q,n-1} \\  
        &        &   \vdots &       &  \vdots &         &   \\ 
        &        &  a'_{n-1,p} &       & a'_{n-1,q} &         &      
\end{bmatrix}
\end{equation}
其中，行数或者列数为p,q的矩阵元发生变化，其他矩阵元不变。
经过计算，矩阵元的变化如下：
\begin{eqnarray}
a'_{pp} &=& c^2 a_{pp} - 2cs a_{pq} + s^2 a_{qq},
\label{jacobi-a'pp}
\\
a'_{qq} &=& s^2 a_{pp} + 2cs a_{pq} + c^2 a_{qq},
\label{jacobi-a'qq}
\\
a'_{pq} &=& (c^2 - s^2)a_{pq} + cs( a_{pp} - a_{qq} ),
\label{jacobi-a'pq}
\\
a'_{ip} &=& c a_{ip} - s a_{iq}, ~~\forall i \neq p,q 
\label{jacobi-a'ip}
\\
a'_{iq} &=& s a_{ip} + c a_{iq}, ~~\forall i \neq p,q.
\label{jacobi-a'iq}
\end{eqnarray}
Jacobi希望让非对角元变为零，所以令$a'_{pq}=0$，即
\begin{equation}
(c^2 - s^2)a_{pq} + cs( a_{pp} - a_{qq} )=0,
\end{equation}
我们就得到了旋转角$\phi$的相关信息
\begin{eqnarray}
\theta = \cot 2 \phi = \frac{ a_{qq} - a_{pp} }{ 2 a_{pq}}.
\end{eqnarray}
如果我们记
\begin{equation}
t = \frac{s}{c} = \tan \phi,
\end{equation}
则有
\begin{equation}
t^2 + 2t\theta - 1 = 0,
\end{equation}
只要满足这个式子，旋转以后就有$a'_{pq}$。
解出来得到
\begin{equation}
t = - \theta \pm \sqrt{\theta^2 +1},
\end{equation}
选取绝对值更小的t，可以使误差更小，
\begin{equation}
t = \left\{
\begin{aligned}
&&- \theta + \sqrt{\theta^2 +1} = \frac{1}{ \theta + \sqrt{\theta^2 +1}}, ~~~ \theta>0, \\
&&- \theta - \sqrt{\theta^2 +1} = - \frac{1}{ -\theta + \sqrt{\theta^2 +1}}, ~~~ \theta<0.
\end{aligned}
\right.
\end{equation}
这两个方程可以用下面的方程概括
\begin{equation}
t = \frac{ sgn(\theta) }{ |\theta| + \sqrt{\theta^2 +1}}.
\end{equation}
在这个约定下，$ |t| < 1$。
所以，我们可以从$a_{pp}, a_{pq}, a_{qq}$计算得到$\theta$，进而用上面的式子得到$t$，得到$t$以后，我们可以进一步计算
\begin{equation}
c = \frac{1}{1+t^2}, ~~ s = ct,
\end{equation}
在上式中，我们假定了$c>0$，综合$|t|<1$，我们实际上约束了$ - \frac{\pi}{4} < \phi < \frac{\pi}{4}$。

出于计算方便和误差控制，我们可以利用$t$，将(\ref{jacobi-a'pp}-\ref{jacobi-a'iq})改写为
\begin{eqnarray}
a'_{pp} &=& a_{pp} - t a_{pq}, \label{jacobi-a'pp-2}, \\
a'_{qq} &=& a_{qq} + t a_{pq}, \label{jacobi-a'qq-2}, \\
a'_{ip} &=& a_{ip} - s (a_{iq} + \frac{s}{1+c} a_{ip}), \\
a'_{iq} &=& a_{iq} + s (a_{ip} - \frac{s}{1+c} a_{iq}).
\end{eqnarray}

到这里，我们理清了旋转变换怎样使$a'_{pq}$等于零。
我们可以进行一系列这样的旋转变换。
有时变为零的非对角元，可能在下一次旋转变换之后，又变为非零。
比如
\begin{equation}
A'' = P^\top_{pr} P^\top_{pq} A P_{pq} P_{pr},
\end{equation}
经过第一次正交变换$P_{pq}$，$a'_{pq}=0$，但是经过第二次正交变换$P_{pr}$，根据(\ref{jacobi-a'ip})，有
\begin{equation}
a''_{pq} = a''_{qp} = ca'_{qp} - sa'_{qq} = -s a'_{qq} \neq 0.
\end{equation}
那么这一系列正交变换会使得A最终变为对角阵吗？
答案是肯定的。

{\bf 任何正交变换前后，矩阵的所有矩阵元平方和不变。}这是正交变换的性质，这很容易证明。
矩阵$N=A^\top A$的对角元为
\begin{equation}
N_{ii} = \sum_j A^\top_{ij} A_{ji} = \sum_j A_{ji} A_{ji},
\end{equation}
所以$N$的迹为
\begin{equation}
Trace(N) = \sum_i N_{ii} = \sum_{ij} A_{ji} A_{ji},
\end{equation}
正是矩阵$A$的所有矩阵元的平方和。
经过正交变换$P$以后，矩阵$A$变为
\begin{equation}
A'=P^\top A P,
\end{equation}
它的所有矩阵元平方和为
\begin{eqnarray}
Trace( A'^\top A) &=& Trace( P^\top A^\top P P^\top A P ) = Trace( P^\top A^\top A P ) 
\nonumber\\
&=& Trace( P^\top N P) = Trace(N) = \sum_{ij} A_{ji} A_{ji}.
\end{eqnarray}
所以正交变换前后，矩阵的所有矩阵元平方和不变。

根据(\ref{jacobi-a'pp}-\ref{jacobi-a'iq})，正交变换$P_{pq}$之后，
\begin{equation}
a'^2_{pp} + a'^2_{qq} = a^2_{pp} + a^2_{pq} + a^2_{qp} + a^2_{qq},
\end{equation}
所以矩阵$A$经过这次变换以后，对角元平方和增加了 $2 a^2_{pq}$。

对角元平方和的上限为$A$的所有矩阵元平方和，而我们可以任意选取$p,q$进行$P_{pq}$变换，所以最终一定可以将$A$对角化。

\iffalse
\begin{eqnarray}
x &=& - a_{pq}, \\
y &=& (a_{qq} - a_{pp})/2, \\
\omega &=& sgn(y)\frac{x}{\sqrt{x^2+y^2}}, \\
s &=& \frac{\omega}{\sqrt{2 + 2\sqrt{1-\omega^2}}} = \frac{ \omega }{ \sqrt{1-\omega} + \sqrt{1+\omega} }, \\
c &=& \sqrt{1-s^2} = \frac{1+\sqrt{1-\omega^2}}{\sqrt{1-\omega} + \sqrt{1+\omega}},
\end{eqnarray}
可以证明，在这样的相似变换后，

所以，$\alpha=(\begin{smallmatrix}
a_{pp} & \\
 & a_{qq} 
\end{smallmatrix})$经过了$\sigma=(\begin{smallmatrix}
c & s \\
-s & c
\end{smallmatrix})$的相似变换，非对角元变为零了。
因为这是一个正交变换，所以矩阵元的平方和不变，这很容易证明：
\begin{eqnarray}
\sum_{i,j} |a'_{ij}|^2 &=& trace( (\sigma^\top \alpha \sigma)^\top (\sigma^\top \alpha \sigma) ) \nonumber\\
&=& trace( \sigma^\top (\alpha^\top \alpha) \sigma ) = trace( \alpha^\top \alpha) \nonumber\\
&=& \sum_{ij} | a_{ij} |^2
\end{eqnarray}
而$\alpha'$的非对角元为零，所以在相似变换以后，对角元的平方和增大了。
而$A$中其它对角元都不变，所以每经过一次这样的相似变换，$A$的对角元平方和就增大了。
正是因为如此，Jacobi方法迭代下去，非对角元的平方和会一直变小，最终收敛为0。
对于行数不为$p,q$的元素，
\begin{eqnarray}
a'_{ip} &=& c a_{ip} + s a_{iq}, \nonumber\\
a'_{iq} &=& -s a_{ip} + c a_{iq},
\end{eqnarray}
对于列数不为$p,q$的元素，
\begin{eqnarray}
a'_{pi} &=& c a_{pi} + s a_{qi}, \nonumber\\
a'_{qi} &=& -s a_{pi} + c a_{qi}.
\end{eqnarray}
\fi
如果把所有相似变换矩阵$P_{pq}$的乘积记作
\begin{equation}
S = P_{p_1 q_1} P_{p_2 q_2} \cdots P_{p_N q_N},
\end{equation}
对角化后的矩阵记为$\Lambda$，则有
\begin{equation}
S^\top A S = \Lambda,
\end{equation}
即
\begin{equation}
AS = S \Lambda,
\end{equation}
也就是说，S的各列是A的特征向量，$\Lambda$的对角元是A的特征值。
所以在程序中，可以一边做正交变换，一边记下它们的总乘积，即同时得到特征值和特征向量。

一般需要$O(n^2)$次旋转变换，Jacobi算法才收敛。
每次旋转变换需要改变2行、2列数字，需要$O(n)$次操作，所以Jacobi算法共需要$O(n^3)$次操作，才会收敛。

下面的代码实现Jacobi对角化。
\begin{lstlisting}
/*
 * jacobi(...): diagonalizes a symmetric real matrix, outputs eigenvalues and eigenkets
 *
 * n is the dimension of the symmetric matrix
 * a is the matrix to be diagonalized, the element Aij=a[i*n+j]
 * z gives the eigenvalues of this matrix A,from small to large
 * v (its column arrays) gives the eigenfunctions corresponding to z. 
 * z's and v's initial values do not matter, they would be assigned in this code
 * e.g. bj=v[j][p](p=0,..,n-1) is the jth(j=0,..,n-1) normalized eigenfuction of z[j].
 */
void jacobi(int n, double *a, double *z, double **v)
{
        int l, i, j, p, q;
        double max, temp, s, c, omega, x, y;
        double g, h;
        for (i = 0; i < n; i++)
        {
                for (j = 0; j < n; j++)
                {
                        v[i][j] = 0;
                }
                v[i][i] = 1.0;         //v is a unit matrix
        }

        for (l = 1; l < max_iteration; l++)
        {
                max = 0;// find the largest off-diagonal matrix element, and its coordinate p,q
                for (i = 0; i < n; i++)
                {
                        for (j = 0; j < n; j++)
                        {
                                if ( i != j && fabs( a[i*n+j] ) >max )
                                {
                                        max = fabs( a[i*n+j] );
                                        p = i;
                                        q = j;
                                }
                        }
                }

                if ( max < 1e-10 )
                        break;

                // construct similarity transformation
                x = -a[p*n+q];
                y = ( a[q*n+q] - a[p*n+p] )/2.0;
                omega = x / sqrt(x*x+y*y);
                if (y < 0) omega = -omega;
                s = 1.0 + sqrt(1.0 - omega*omega);
                s = omega / sqrt(2.0*s);
                c = sqrt(1.0 - s*s);
                temp = a[p*n+p];
                a[p*n+p] = temp * c * c + a[q*n+q] *s *s + a[p*n+q] *omega;
                a[q*n+q] = temp * s * s + a[q*n+q] *c *c - a[p*n+q] *omega;
                a[p*n+q] = 0;
                a[q*n+p] = 0;

                for (j = 0; j < n; j++)
                {
                        if (j != p && j != q)
                        {
                                temp = a[p*n+j];
                                a[p*n+j] = temp * c + a[q*n+j] * s;
                                a[q*n+j] = -temp * s + a[q*n+j] * c;
                        }
                }
                for (i = 0; i < n; i++)
                {
                        if (i != p && i != q)
                        {
                                temp = a[i*n+p];
                                a[i*n+p] =  temp * c + a[i*n+q] * s;
                                a[i*n+q] = - temp * s + a[i*n+q] * c;
                        }
                }
                for (i = 0; i < n; i++)// v = S_1 S_2 ..., its columns are eigenvectors of matrix *a
                {
                        temp = v[i][p];
                        v[i][p] =  temp * c + v[i][q] * s;
                        v[i][q] = - temp * s + v[i][q] * c;
                }
        }//l
        if( l== max_iteration ){
                cout<<" Jacobi: after "<<max_iteration<<" iterations, didn't converge...\n";
                exit(1);
        }
        for (i = 0; i < n; i++)
                z[i] = a[i*n+i];
        for(p=0;p<n-1;p++)// sort the eigenvalues and eigenfunctions, so that eigenvalues are arranged from small to large
        {
                temp=z[p];
                l=p;
                for(i=p+1;i<n;i++){// bubble sorting
                        if(g>z[i]){// find the smallest eigenvalue in z[p], ..., z[n-1]
                                temp=z[i];
                                l=i;
                        }
                }
                if(l!=p)
                {
                        temp=z[p];// put the smallest eigenvalue in {z[p],...,z[n-1]} in z[p]
                        z[p]=z[l];
                        z[l]=temp;
                        for(j=0;j<n;j++)// exchange v[][p] with v[][l]
                                // this can be optimized, we can store eigenwave in a row, rather than in a column
                        {
                                temp = v[j][p];
                                v[j][p]=v[j][l];
                                v[j][l]=temp;
                        }
                }
        }                      //eigvalues from small to large
        for(i=0;i<n;i++)
        {
                temp=0;
                for(j=0;j<n;j++)
                {
                        temp+=v[j][i]*v[j][i];
                }
                for(j=0;j<n;j++)
                {
                        v[j][i] /= sqrt(temp);  //normalization of wavefunctions
                }
        }
}
\end{lstlisting}

\section{将实对称阵三对角化}
在上一章中，我们介绍了追赶法解三对角线性方程组，只需要$O(n)$次操作。
对于一般的线性方程组，用高斯消元法、或者LU分解，需要$O(n^3)$次操作。

所以，如果能把矩阵转化为三对角阵，问题一般就好解决了，实对称阵对角化，也是如此。
下面我们介绍两种三对角化的方法；一种是Givens变换，用一系列特殊设计的转动矩阵，实现三对角化；另一种是Householder变换，用一系列镜面反射变换，实现目标。

\subsection{*Givens变换}
Givens变换完全基于转动变换，我们把上一节关于转动变换的几个公式再放到这里，方便您的阅读。

关于第$p,q$行、第$p,q$列的转动变换矩阵为
\begin{equation}
P_{pq} = \begin{bmatrix}
 1 \\
  & \ddots \\
  &     & c     &     & -s \\
  &     &       & \ddots \\
  &     & s    &         & c \\
  &     &       &        &       & \ddots \\
  &     &       &        &       &       & 1
\end{bmatrix}
\end{equation}
在旋转变换以后，矩阵A变为
\begin{equation}
A' = P^\top_{pq} A P_{pq}
= \begin{bmatrix}
        &        &  a'_{0p} &       & a'_{0q} &         &   \\
        &        &   \vdots &       &  \vdots &         &   \\
a'_{p0} & \cdots & a'_{pp}  & \cdots& a'_{pq} &\cdots   & a'_{p,n-1} \\
        &       & \vdots    &       & \vdots  &     \\
a'_{q0} & \cdots & a'_{qp}  & \cdots& a'_{qq} &\cdots   & a'_{q,n-1} \\  
        &        &   \vdots &       &  \vdots &         &   \\ 
        &        &  a'_{n-1,p} &       & a'_{n-1,q} &         &      
\end{bmatrix}
\end{equation}
也就是说，$A$中，行数与列数都不是$p,q$的矩阵元不变，行数或者列数沾了$p,q$的，就要改变。

在Jacobi那一节中，我们通过调节$c,s$，使得$a‘_{pq}=0$。
这里我们不那么做，我们设计$c,s$，使得$P_{12}$变换之后，$a'_{02}=0$，相应地，也有$a'_{20}=0$，因为$A'$是对称的。
然后，我们设计一个$P_{13}$，使得变换后$a'_{03}=a'_{30}=0$，注意，因为$a'_{02}$不沾$1,3$，所以会保持为零。
如此这般，我们依次设计$P_{14},P_{15},\cdots,P_{1,n-1}$，使得$a'_{04}, \cdots, a'_{0,n-1}$都等于零，经过这$n-1$次转动变换，将第一行、第一列除了$a'_{00}, a'_{01}=a'{10} \neq 0$，其他元素都变为零，所以就朝三对角化靠近了一步。

然后，我们设计一个$P_{23}$，使得$a'_{13}=a'_{31}=0$，那么，$a_{02}, a_{03}$会受影响，重新变为非零吗？
答案是不会。

将矩阵A右边乘一个矩阵，实则是对A各列进行线性变换，即乘了以后，各列是A原来各列的线性叠加。
而将矩阵A左边乘一个矩阵，是对A各行进行线性变换。
所以，做变换$P_{23}$，是对A进行第2,3列的线性重组，然后对第2,3行进行线性重组。
变换以后，$a'_{02},a'_{03}$是\footnote{这里的$a$和$a'$只是表示一次旋转变换前后的矩阵，所以不同变换对应的$a,a'$是不一样的，如果引起混乱，抱歉。}$a_{02}, a_{03}$的线性叠加，而这次$P_{23}$变换之前，$a_{02}=a_{03}=0$，所以仍有$a'_{02}=a'_{03}=0$。
所以，我们依次进行$P_{23}, P_{24}, \cdots, P_{2,n-1}$，将$a_{13}, \cdots, a_{1,n-1}$都变为零，这些操作之后，第一行后面$n-2$个元素仍然保持为零。

照这个道理进行下去，所有转动变换的乘积构成Givens变换，
\begin{equation}
S = P_{12}\cdots P_{1,n-1} P_{23} \cdots P_{2,n-1} \cdots P_{n-2,n-1}.
\end{equation}
在这些变换下，原来的实对称阵会变为三对角阵。

每次旋转变换$P_{pq}$，需要改变$A$中第$p,q$行、第$p,q$列，共进行$O(n)$次操作。
Givens变换需要$O(n^2)$次

\subsection{Householder变换}
Householder变换是一种正交变换，通过它可以将对称阵转化为三对角矩阵，三对角矩阵就容易处理了。

Householder变换由一系列正交相似变换构成,
\begin{eqnarray}
S &=& S_1 S_2 \cdots S_{n-2}, \\
S_1 S^\top_1 &=& S_2 S^\top_2 = \cdots = 1, \\
S^\top A S &=& 
\begin{pmatrix}
a_{11} & e_2 &    &  & \\
e_2 & a'_{22} & e_2 &  & \\
    & \ddots & \ddots & \ddots \\
    &     &         & a^{(n-2)}_{n-1,n-1} & e_{n-1} \\
    &     &         & e_{n-1}             & a^{(n-2)}_{n,n}
\end{pmatrix}
\end{eqnarray}
上面的$a'_{22}$经过了一次迭代（一次相似变换），$a^{(n-2)}_{n-1,n-1}, a^{(n-2)}_{n,n}$经过了$n-2$次迭代。

$S_1$这样构造：
\begin{eqnarray}
S_1 = \begin{pmatrix}
1 & \vec{0}^\top \\
\vec{0} & P
\end{pmatrix},
\end{eqnarray}
其中$\vec{0}$表示n-1个0构成的列矢量，$\vec{0}^\top$是它的转置，$P$满足$P^\top = P, P^2 = 1$。
这样构造：
\begin{equation}
P = 1 - 2 \vec{u} \vec{u}^\top,
\end{equation}
其中$\vec{u}$是归一化的矢量，那么P的这个形式就可以保证$P^\top = P, P^2 = 1$。

为了方便叙述，我们把A分为四块矩阵
\begin{equation}
A = \begin{pmatrix}
a_{11} & \vec{v}^\top \\
\vec{v} & B
\end{pmatrix},
\end{equation}
其中，$\vec{v} = \{ a_{21}, \cdots, a_{n1} \}^\top$。

经过第一次相似变换以后，$A$变为
\begin{equation}
S^\top_1 A S_1 =
\begin{pmatrix}
a_{11} & (P \vec{v})^\top \\
P\vec{v} & B'
\end{pmatrix},
\end{equation}
其中$B' = PAP$。

所以我们想，如果$P\vec{v}$只有第一个元素不为零，其它都是零，我们就离目标（三对角阵）接近了一步。

$(P \vec{v})$的模平方为
\begin{equation}
(P \vec{v})^\top (P\vec{v}) = |\vec{v}|^2,
\end{equation}
所以，我们得命令
\begin{equation}
P \vec{v} = \pm v \vec{e}_2,
\end{equation}
其中v是$\vec{v}$的模，$\vec{e}_2$是$\vec{v}$的线性空间第一个单位矢量。

这样，我们就把$\vec{u}$定下来了，
\begin{equation}
P\vec{v} = (1- 2 \vec{u} \vec{u}^\top) \vec{v}
= \vec{v} - 2 \vec{u} (\vec{u}^\top \vec{v}) 
= \pm v \vec{e}_2,
\end{equation}
最后一个等号两边同时与$\vec{v}$做点乘，得到
\begin{equation}
2(\vec{u}^\top \vec{v} )^2 = v^2 \mp a_{21} v.
\end{equation}
所以可以定义
\begin{equation}
\vec{u} = \frac{ \vec{v} \mp v \vec{e}_2 }{ 2 \vec{u}^\top \vec{v} } 
= \frac{\vec{v} \mp v \vec{e}_2}{\sqrt{ 2v^2 \mp 2a_{21} v  }},
\end{equation}
其中$\mp$对应两种选择，都可以完成目标，但我们一般选$v^2 \mp a_{21} v$更大的那一种，防止误差太大。

\begin{shaded*}
这里我们可以做一个观察，如果把$\vec{u}$定为
\begin{equation}
\vec{u} = \frac{\vec{v} \mp v \vec{e}_k}{\sqrt{ 2v^2 \mp 2a_{k1} v  }}, ~~~ k=2,3,...,n
\end{equation}
那么，$P=1-2\vec{u}\vec{u}^\top$作用到$\vec{v} = \{ a_{21}, a_{31}, \cdots, a_{n1} \}^\top$上，会得到
\begin{equation}
P\vec{v} = \vec{v} - 2 \vec{u} (\vec{u}^\top \vec{v}) = \pm v \vec{e}_k.
\end{equation}
也就是，$S^\top_1 A S_1$第一列只有第1个和第k个元素不为零，第一行只有第1个和第k个元素不为零。
所以，这个技巧可以针对$k=2,...,n$的，使得$\vec{v}$中的相应的分量不为零，其余变为零。
这样对一个矢量$\vec{v}$做变换，原本就叫做Householder变换，它实际上是将矢量$\vec{v}$相对一个超平面做镜像反射变换，这个超平面的法向矢量为$\vec{u}$。

我们可以再做一个观察，如果把$S_1$定为
\begin{eqnarray}
S_1 = \begin{pmatrix}
P & \vec{0}^\top \\
\vec{0} & 1
\end{pmatrix},
\end{eqnarray}
矩阵A划分为
\begin{equation}
A = \begin{pmatrix}
B & \vec{v} \\
\vec{v}^\top & a_{nn}
\end{pmatrix},
\end{equation}
会发现
\begin{equation}
\begin{pmatrix}
A' & P \vec{v}\\
(P\vec{v})^\top & a_{nn}
\end{pmatrix},
\end{equation}
其中$\vec{v} = \{ a_{1n}, a_{2n}, \cdots, a_{n-1,n} \}$。
然后，我们可以定义
\begin{equation}
\vec{u} = \frac{\vec{v} \mp v \vec{e}_{n-1}}{\sqrt{ 2v^2 \mp 2a_{n-1,1} v  }}, ~~~ k=2,3,...,n
\end{equation}
这样，$S^\top_1 A S_1$的最后一行、最后一列就只有第$n-1,n$个元素不为零，其它都是零。

所以，我们这一节介绍的Householder变换，也可以从最后一行（列）开始，逐渐向三对角形式靠近。
\end{shaded*}

把定下的$\vec{u}$带入P，就得到了$S_1$，然后进行相似变换，就完成了第一步。
但是直接做两次矩阵乘法，太费时间了，我们可以如下简化计算
\begin{eqnarray}
\vec{p} &=& 2 B \vec{u}, \\
BP &=& B(1-2\vec{u}\vec{u}^\top) = B - \vec{p} \vec{u}^\top, \\
PBP &=& (1-2\vec{u}\vec{u}^\top)(B - \vec{p} \vec{u}^\top) \nonumber\\
    &=& B - \vec{u}\vec{p}^\top - \vec{p}\vec{u}^\top + 2 (\vec{u}^\top \vec{p}) \vec{u}\vec{u}^\top.
\end{eqnarray}
如果定义
\begin{eqnarray}
K &\equiv& \vec{u}^\top \vec{p}, \nonumber\\
\vec{q} &=& \vec{p} - K \vec{u},
\end{eqnarray}
则有
\begin{eqnarray}
PBP = B - \vec{u}\vec{q}^\top - \vec{q}\vec{u}^\top.
\end{eqnarray}

说到这里，我们解释了$S_1$的构造思路，并且展示了效果以及计算方法。
$S_2$则只作用在$B'$上，将$B'$第一行、第一列第$3,\cdots,n-1$个元素变为零。
如此继续下去，通过$n-2$次相似变换，最终将$A$转换为三对角形式。

下面的代码用n-2次householder变换，将一个实对称阵三对角化，并记下所有相似变换矩阵的乘积
\begin{equation}
S = S_0 S_1 \cdots S_{n-3},
\end{equation}
S使得
\begin{equation}
S^\top A S = (tridiag).
\end{equation}
(tridiag)表示三对角形式的矩阵。
在下面的代码中，S被写入a，覆盖掉a原来的数据（矩阵A）。

\begin{lstlisting}
/* householder(n,a,d,e):        tridiagonalize a real symmetric matrix by householder transformation
 * n:           dimension of the matrix
 * a[n*n]       the original matrix, will be overwritten in this function, end up to be the similarity transformation matrix S = S_0 ... S_{n-3}, so that S^top A S = tridiag
 * d[n]         the diagonal elements after tridiagonalization
 * e[n]         the sub-diagonal elements after tridiagonalization, e[0] will be set to 0.
 *
 * After tridiagonalization:
 * d0 e1
 * e1 d1 e2
 *    e2 d3 e3
 *       .......
 *          .......
 *             e_{n-1} d_{n-1}
 */
void householder(int n, double *a, double *d, double *e){
        int i,j,k,m;
        double v, denorm, K, y;
        double *u, *p, *q;
        u = new double [n];
        p = new double [n];
        q = new double [n];

        e[0]=0;

        for(i=0;i<n-2;i++){// tridiagonalize the 0, 1, ..., n-2 th column & row, succesively
                // ***************** calculate \vec{v}, \vec{u}
                v=0;
                for(j=i+1;j<n;j++){
                        u[j] = a[i*n+j];// \vec{u} = \vec{v} temporarily
                        v += a[i*n+j] * a[i*n+j];
                }
                v = sqrt(v);// sqrt( a^2_{i,i+1} + ... + a^2_{i,n} )
                if( fabs(v) < 1E-14){// skip one transformation
                        d[i]=a[i*n+i]; e[i+1] = 0;
                        continue;
                }

                denorm = 2*v*v;// denorm = sqrt( 2v^2 -+ 2a_{i+1,i} v )
                if( a[i*n+i+1] > 0 ){
                       denorm += 2* a[i*n+i+1] * v;
                       a[i*n+i+1] = - v;// P\vec{v} = +- v \vec{e}_i
                       a[(i+1)*n+i] = -v;
                }
                else{
                        denorm += - 2* a[(i+1)*n+i] *v;
                        a[i*n+i+1] = v;
                        a[(i+1)*n+i] = v;
                }
                d[i]=a[i*n+i]; e[i+1] = a[i*n+i+1];// record d[i] and e[i+1]
                for(j=i+2;j<n;j++){// tridiagonalization of the ith column & row
                        a[i*n+j] = 0;
                        a[j*n+i] = 0;
                }
                denorm = sqrt(denorm);// denorm = sqrt( 2v^2 -+ 2a_{i+1,i} v )

                for(j=i+1;j<n;j++){// \vec{u} = \vec{v}/denorm -+ v/denorm \vec{e}_i
                        u[j] /= denorm;
                }
                if( u[i+1] > 0 ) u[i+1] += v/denorm;
                else u[i+1] -= v/denorm;

                // ***************** calculate K
                K=0;//calculate K = \vec{u}^\top \vec{p}
                for(j=i+1;j<n;j++){// determine p
                        y=0;
                        for(k=i+1;k<n;k++) y += a[j*n+k] * u[k];
                        p[j] = 2*y;// \vec{p} = 2 B \vec{u}
                        K += p[j] * u[j];
                }

                // ***************** calculate \vec{q}
                for(j=i+1;j<n;j++){// determine \vec{q} = \vec{p} - K \vec{u}
                        q[j] = p[j] - K*u[j];
                }

                // ***************** calculate PBP
                for(j=i+1;j<n;j++){// PBP = B - \vec{u}\vec{q}^\top - \vec{q}\vec{u}^\top
                        for(k=i+1;k<n;k++){
                                a[j*n+k] -= u[j]*q[k] + q[j]*u[k];
                        }
                }

                for(j=i+1;j<n;j++) a[ i*n + j ] = u[j];//restore u, to calculate P0P1P2P3...P_{n-3}
                /* for test only
                cout<<"i="<<i<<"\t A = "<<endl;
                for(j=0;j<n;j++){
                        for(k=0;k<n;k++){
                                cout<<a[j*n+k]<<",";
                        }
                        cout<<endl;
                }
                */
        }
        d[n-2]=a[(n-2)*n+n-2]; d[n-1] = a[(n-1)*n+n-1];
        e[n-1]=a[(n-1)*n+n-2];
        // ************** write P_0 P_1 ... P_{n-3} into matrix *a
        for(j=n-2;j<n;j++){// set the right-lowest 2x2 matrix as unit matrix
                for(k=n-2;k<n;k++){
                        if(k==j)a[j*n+k]=1;
                        else a[j*n+k]=0;
                }
                cout<<endl;
        }
        for(i=n-3;i>=0;i--){// calculate S=S_0 S1...S_{n-3} and stores it in a
                for(j=i+1;j<n;j++){// \vec{u}^\top B is stored in A[i][i+1]...A[i][n-1]
                        y=0;// y= 2\sum_k u_k a_{kj}
                        for(k=i+1;k<n;k++)
                                y+= a[i*n+k]*a[k*n+j];
                        y*=2;
                        for(k=i+1;k<n;k++){// a_{ij} = a_{ij} - 2 u_i \sum_k u_k a_{kj}
                                a[k*n + j] -= a[i*n+k]*y;
                        }
                }
                for(j=i+1;j<n;j++){
                        a[i*n+j]=0;// delete info of \vec{u}
                        a[j*n+i]=0;
                }
                a[i*n+i]=1;
        }

        delete [] u,p,q;
}
\end{lstlisting}
下面是一个主函数，测试上面的功能函数。
定义一个随机的对称矩阵a\_house，将方阵a\_house三对角化$S^\top (a\_{house}) S = (tridiag)$，然后再计算$S(tridiag)S^\top$，看看是否等于a\_house，如果等于，则说明程序的计算是正确的。
\begin{lstlisting}
int main(){

        srand(time(NULL));

        int n=10,i,j,k;
        double tempdouble,y;
        double *a_house= new double [n*n];
        cout<<"The original matrix: "<<endl;
        for(i=0;i<n;i++){
                for(j=i;j<n;j++){
                        tempdouble = (double)rand()/RAND_MAX;
                        a_house[i*n+j] = tempdouble;
                        a_house[j*n+i] = tempdouble;
                }
                for(j=0;j<n;j++){
                        cout<<a_house[i*n+j]<<",";
                }
                cout<<endl;
        }
        cout<<"ready? Go!"<<endl;

        double *d= new double [n];
        double *e= new double [n];
        householder(n,a_house,d,e);

        double *a_mid = new double [n*n];
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        y=a_house[i*n+j-1]*e[j] + a_house[i*n+j]*d[j];
                        if(j<n-1) y += a_house[i*n+j+1]*e[j+1];
                        a_mid[i*n+j] = y;
                }
        }
        cout<<"a_house * tridiag * a_house^top: "<<endl;
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        y=0;
                        for(k=0;k<n;k++)
                                y += a_mid[i*n+k] * a_house[j*n+k];
                        cout<<y<<",";
                }
                cout<<endl;
        }

        delete [] a_house, z, d, e;

        return 0;
}
\end{lstlisting}

\section{三对角矩阵对角化}

\paragraph{QR分解}
每个实对称阵$A$，都可以表示成
\begin{equation}
A = QR,
\end{equation}
其中，$Q$是一个正交矩阵，$R$是一个上三角矩阵。
前人证明了\footnote{我还没有学习这个证明，也许过两年会有时间看明白，但现在不妨碍我们实践。}，若反复进行正交变换
\begin{equation}
A' = RQ = Q^\top A Q,
\end{equation}
最终$A'$会收敛为对角阵，本征值按绝对值从小到大排列，如果有个特征值$\lambda_i$有$p_i$重简并，则会出现一个$p_i$维子矩阵出现在对角线上，这个子矩阵的特征值均为$\lambda_i$。

根据$A$找到$Q$，是可以做到的，比如，我们可以用上一节介绍的技巧，构造一系列Householder变换，左乘到A上，使得A的下三角部分归零，
\begin{equation}
P_{n-2} \cdots P_0 A = R,
\end{equation}
定义
\begin{equation}
Q^\top =  P_{n-2} \cdots P_0,
\end{equation}
则有
\begin{equation}
A = QR,
\end{equation}
即完成QR分解。

\paragraph{QL分解}
类似地，我们也可以构造Householder变换$P_{n-1}$，使得A的最右列除最后一个元素外，其余元素归零，类似地构造$P_{n-2}, \cdots, P_{1}$，使得
\begin{equation}
P_1 \cdots P_{n-1} A = L,
\end{equation}
其中$L$表示下三角矩阵。我们定义
\begin{equation}
Q^\top =  P_{1} \cdots P_{n-1},
\end{equation}
则有
\begin{equation}
A = Q L,
\end{equation}
完成QL分解。
若反复做相似变换
\begin{equation}
A' = Q^\top A Q,
\end{equation}
最终也会收敛为对角阵，或包括对角子矩阵块（若有简并）的对角阵。

古人证明了，QR或QL变换，不改变矩阵的三对角性，即每次变换以后，A'仍然是三对角阵。

但是，每一次正交变换，需要$O(n)$次householder变换，而一次householder变换需要$O(n^2)$次操作，所以一次这样的正交变换，就需要$O(n^3)$次操作，还不知道要几次正交变换才能收敛。
所以，这个办法对于一般的实对称阵的对角化，是非常不人道的。
但是，对于三对角阵，每次householder的变换矩阵只有2x2个元素与单位阵不一样，所以$O(n)$次householder变换，所需要的总操作数是$O(n)$。
所以，QR或者QL变换用来对角化三对角矩阵，是有效的。

\paragraph{所有本征值进行位移}
根据前人的说法\cite{NumericalRecipes}，若$|\lambda_j| > |\lambda_i|$，s次迭代以后，
\begin{equation}
a^{(s)}_{ij} \sim ( \frac{\lambda_i}{\lambda_j} )^s,
\end{equation}
所以$\lambda_i / \lambda_j$越小，收敛得越快。
如果让A减去单位阵的常数，那么A的本征值会全部减去这个常数，即所有本征值进行一个平移
如果平移量很接近$\lambda_i$，平移以后，$\lambda_i$变成一个很小的量，$\lambda_i / \lambda_j$就是一个小量，QR或者QL变换收敛的就快。
所以这个平移的技巧会加速收敛。

也就是说
\begin{equation}
A = QL + k_s E,
\end{equation}
其中$k_s$为位移量，那么
\begin{equation}
A' = LQ + k_s E = Q^\top A Q,
\end{equation}
相比不用位移，得到的$Q$是不同的。

根据古人的说法，可以这么做
\begin{shaded*}
\begin{itemize}
\item [1] 求解三对角矩阵左上角2x2的小矩阵的本征值$k_s$，选取距离$a_{00}$更接近的那个，按上面的公式进行位移，然后进行一系列QL变换，$a_{01}$会很快收敛。
\item [2] 然后求解$(\begin{smallmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22}
\end{smallmatrix})$的本征值$k_s$，选取距离$a_{11}$更接近的那个，然后位移，进行一系列QL变换，$a_{12}$会很快收敛。
\item[3] 然后求解$(\begin{smallmatrix}
a_{22} & a_{23} \\
a_{32} & a_{33}
\end{smallmatrix})$的本征值$k_s$。。。
\end{itemize}
\end{shaded*}
若求$(\begin{smallmatrix}
a_{pp} & a_{pq} \\
a_{qp} & a_{qq}
\end{smallmatrix})$的本征值，会得到
\begin{equation}
x = \frac{ a_{pp} + a{qq} }{2} \pm \sqrt{(\frac{a_{pp} - a{qq}}{2})^2 + a^2_{qq}},
\end{equation}
选取其中离$a_{pp}$更近的根，则为
\begin{eqnarray}
k_s &=& \frac{ a_{pp} + a{qq} }{2} - sgn(a_{qq} - a_{pp}) \sqrt{(\frac{a_{pp} - a{qq}}{2})^2 + a^2_{qq}} \nonumber \\
&=& a_{pp} - sgn(\theta) \frac{a_{pq} }{ |\theta| + \sqrt{\theta^2 +1} },
\end{eqnarray}
其中$\theta = (a_{qq}-a_{pp})/(2a_{pq})$，如同在Jacobi对角化一节定义的一样。

根据古人的经验，这样引入位移进行QL变换，这个三对角矩阵会很快收敛为对角阵，QL变换的次数为$O(n)$。

\paragraph{QL分解的自由度}
若有
\begin{equation}
A' = Q^\top A Q,
\end{equation}
而$A'$是一个三对角阵，可以记作
\begin{equation}
A' = \begin{bmatrix}
d_1 & e_{2} \\
e_2 & d_2 & d_3\\
    &\ddots & \ddots & \ddots \\
    &		& e_{n-2} & d_{n-1} & e_{n-1} \\
    &		&		& e_{n-1} & d_n 
\end{bmatrix}
\end{equation}
然而，
\begin{equation}
A'Q^\top = Q^\top A,
\end{equation}
把Q分为$Q=(\vec{q}_1, \cdots, \vec{q}_n)$，
所有上式等号两边最后一行为
\begin{equation}
e_{n-1} \vec{q}^\top_{n-1} + d_{n} \vec{q}^\top_n = \vec{q}^\top_n A,
\end{equation}
{\bf 如果$\vec{q}_n$已经确定了}，那么根据上式，可以唯一地确定
\begin{equation}
d_n = \vec{q}^\top_n A \vec{q}_n,
\end{equation}
以及
\begin{equation}
e_{n-1} \vec{q}^\top_{n-1} = \vec{q}^\top_n A - d_{n} \vec{q}^\top_n,
\end{equation}
$\vec{q}_{n-1}$是归一化的，所以唯一的变数是：$e_{n-1}$的$\pm$符号，相应的$\vec{q}_{n-1}$是$\mp$符号。
$A'Q^\top = Q^\top A$的倒数第二行为
\begin{equation}
e_{n-2} \vec{q}^\top_{n-2} + d_{n-1} \vec{q}^\top_{n-1} + e_{n-1} \vec{q}^\top_{n} = \vec{q}^\top_{n-1} A,
\end{equation}
得到
\begin{equation}
d_{n-1} = \vec{q}^\top_{n-1} A \vec{q}_{n-1},
\end{equation}
不管$\vec{q}_{n-1}$的符号如何，上式的值都是唯一确定的。
另有
\begin{equation}
e_{n-1} (e_{n-2} \vec{q}_{n-2}^\top) = (e_{n-1} \vec{q}_{n-1})^\top A - d_{n-1} (e_{n-1}\vec{q}_{n-1}^\top) - e^2_{n-1} \vec{q}^\top_n,
\end{equation}
所以$e_{n-1} e_{n-2} \vec{q}_{n-2}^\top$是唯一确定的。

如此类推，会发现，$\vec{q}_n$确定以后，$Q^\top$基本上就确定了，除了各行的符号有不确定度，而这些不确定度对应着不同但是都成立的QL分解。

\paragraph{构造QL分解}
可以使用Householder变换，但是很复杂，不如用旋转变换。
因为$Q^\top A = L$，所以可以如下构造QL分解
\begin{shaded*}
\begin{itemize}
\item [1] 先构造一个旋转变换$P_{n-1,n}$，左乘到A以后，使得$a_{n-1,n}$变为零。
\item [2] 构造一个旋转变换$P_{n-2,n-1}$，左乘到$P_{n-1,n}A$，使得$a_{n-2,n-1}$变为零。
\item [3] 。。。
\item [4] 构造一个旋转变换$P_{1,2}$，使得$a_{1,2}$变为零。
\end{itemize}
\end{shaded*}
这样会得到
\begin{eqnarray}
P_{1,2}\cdots P_{n-1,n} A &=& L, \\
A' &=& P_{1,2} \cdots P_{n-1,n} A P^\top_{n-1,n} \cdots P^\top_{1,2}.
\end{eqnarray}

也可以用另一种方法构造。
前面我们说明了，如果$Q^\top$的最后一行定下来，那么它基本上就定下来了，除了各行的正负号具有任意性。
我们注意到$P_{1,2} \cdots P_{n-1,n}$的最后一行，实际上就是$P_{n-1,n}$的最后一行。
这意味着，我们可以用其他办法构造$Q^\top$，只要最后一行与$P_{n-1,n}$的最后一行相同，我们得到的$Q$就与$P_{1,2} \cdots P_{n-1,n}$基本一致，最多各行差一些正负号，而且也构成QL分解。
\begin{shaded*}
\begin{itemize}
\item[1] 构造$P_{n-1,n}$旋转变换，使得$a_{n-1,n}$变为零。这会使得$a_{n-2,n}$变为非零，其他保持三对角形式。
\item[2] 构造$P_{n-2,n-1}$旋转变换，使得$a_{n-2,n}$变为零。这会使得$a_{n-3,n-1}$变为非零，其他保持三对角形式。
\item[3]。。。
\item[n-1] 构造$P_{1,2}$旋转变换，使得$a_{1,3}$变为零。回复三对角化。
\end{itemize}
\end{shaded*}
用这个算法，会得到
\begin{equation}
A' = P^\top_{1,2} \cdots P^\top_{n-1,n} A P_{n-1,n} \cdots P_{1,2},
\end{equation}
我们只需依次计算
\begin{eqnarray}
&&P^\top_{n-1,n} A P_{n-1,n} \nonumber\\
&&P^\top_{n-2,n-1}P^\top_{n-1,n} A P_{n-1,n}P_{n-2,n-1}, \nonumber\\
&&\cdots \nonumber
\end{eqnarray}

可以加入位移，写成代码，如下
\lstset{basicstyle=\scriptsize}
\begin{lstlisting}
/*
 * tqls: diagonalize a tridiagonal matrix in QL algorithm with shifting
 *
 * int n:               dimension
 * double *d            diagonal elements of the tridiagonal matrix, end up to be the eigen values
 * double *e            subdiagonal elements of the tridiagonal matrix
 * double **v           input as [n][n] matrix, whatever value it has, end up to be the eigen kets
 *
 * convention: e[0]=0
 *
 * d[0] e[1]
 * e[1] d[1]
 * ...
 *                      e[n-1]
 *              e[n-1]  d[n-1]
 */
void tqls(int n, double *d, double *e, double **v){

        int m,l,i,iter,j;
        double c,s,r,theta,ks,app,aqq,apq,butterfly;

        // set v as elementary matrix
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        if(i==j) v[i][j]=1;
                        else v[i][j]=0;
                }
        }

        for(l=0;l<n-1;l++){
                iter=0;//迭代使得e[l+1]足够小，需要的步数
                do{
                        for(m=l;m<n-1;m++)//e[m+1]是否已经足够小
                                if( abs(e[m+1]) <= 1E-9*(abs(d[m])+abs(d[m+1])))
                                        break;
                        //执行完上面的代码块，若e[l+1]足够小，则m=l；若e[l+p]足够小，则m=l+p-1;若都没有足够小，则m=n-1
                        //下面的代码处理 l, ..., m 子矩阵块
                        if(m!=l){//e[l+1]还没有足够小
                                if(iter++ == 30){
                                        cout<<"Too many iterations in tqli, when trying to decrease e["<<l<<"] l=0,...,n-1 \n";
                                }
                                theta = (d[l+1]-d[l])/(2.0*e[l+1]);
                                ks = - e[l+1]/(fabs(theta) + sqrt(theta*theta+1));
                                if( theta < 0 ) ks *= -1;
                                ks += d[l];// ks is set up, shifting

                                for(i=l;i<=m;i++) d[i] -= ks;

                                app = d[m-1]; aqq = d[m]; apq = e[m];
                                r = sqrt( apq * apq + aqq * aqq);// c * apq - s *aqq = 0
                                c = aqq / r;
                                s = apq / r;
                                // these codes heavily rely on formulas about Jacobi transformation on my Lecture notes, or on Numerical Recipes 
                                d[m-1] = c*c* app - 2*c*s*apq + s*s*aqq;
                                d[m] = s*s*app + 2*c*s*apq + c*c*aqq;
                                e[m] = (c*c-s*s)*apq + c*s*(app-aqq);
                                butterfly = s* e[m-1];// a'_{m-2,m} = a'_{m,m-2} = -s a_{m,m-1}
                                e[m-1] = c* e[m-1];// a'_{m-2,m-1} = c*a_{m-2,m-1}

                                for(j=0;j<n;j++){
                                        app = v[j][m-1];
                                        aqq = v[j][m];
                                        v[j][m-1] = c*app - s*aqq;
                                        v[j][m] = s*app + c*aqq;
                                }

                                for(i=m-1;i>=l+1;i--){
                                        r = sqrt( butterfly*butterfly + e[i+1]*e[i+1]);
                                        s = butterfly / r;
                                        c = e[i+1] / r;
                                        //butterfly = a_{i+1,i-1}, 
                                        //butterfly' = a'_{i+1,i-1} = c butterfly - s a_{i+1,i} = c butterfly - s e[i+1] = 0

                                        app = d[i-1]; aqq=d[i]; apq = e[i];
                                        d[i-1] = c*c* app - 2*c*s*apq + s*s*aqq;
                                        d[i] = s*s*app + 2*c*s*apq + c*c*aqq;
                                        e[i] = (c*c-s*s)*apq + c*s*(app-aqq);
                                        e[i+1] = s*butterfly + c*e[i+1];
                                        butterfly = s*e[i-1];
                                        e[i-1] = c*e[i-1];

                                        for(j=0;j<n;j++){
                                                app = v[j][i-1];
                                                aqq = v[j][i];
                                                v[j][i-1] = c*app - s*aqq;
                                                v[j][i] = s*app + c*aqq;
                                        }

                                }
                                for(i=l;i<=m;i++) d[i] += ks;
                        }
                }while( m!=l );
        }
}
\end{lstlisting}

\section{Lanczos方法}
Lanczos在大规模稀疏矩阵对角化问题中非常有效。
在原子核壳模型中，Lanczos方法使我们可以处理维数为$10^{10}\sim10^{11}$的大型稀疏矩阵，从而进行一些原子核从头计算。

\subsection{Krylov子空间}
若对实对称阵A，定义\cite{Golub2013Matrix} Rayleigh商
\begin{equation}
r_A( \vec{x} ) = \frac{ \vec{x}^\top A \vec{x} }{ \vec{x}^\top \vec{x} }, ~~\vec{x} \neq 0.
\label{Rayleigh1}
\end{equation}
那么，$\vec{x}$可以表示成A的本征矢量的线性叠加
\begin{eqnarray}
\vec{x} &=& \sum^{n}_{i=1} c_i \vec{a}_i, \\
A \vec{a}_i &=& \lambda_i \vec{a}_i, i=1,...,n, 
\end{eqnarray}
则有
\begin{equation}
r_A( \vec{x} ) = \frac{ \sum_{i,j} c_i c_j \vec{a}^\top_i A \vec{a}_j }{ \sum_i c^2_i }
= \frac{ \sum_i c^2_i \lambda^2_i }{ \sum_i c^2_i },
\label{Rayleigh2}
\end{equation}
所以有
\begin{equation}
\lambda_{min} (A) \leq r_A (\vec{x}) \leq \lambda_{max} (A),
\label{eqn:Reighleigh-2}
\end{equation}
其中$\lambda_{min} (A)，\lambda_{max}(A)$表示A的最小、最大本征值。

也就是说，任意非零的矢量$\vec{x}$处，$A$的Reileigh商都不超出闭区间$[\lambda_{min}(A), \lambda_{max}(A)]$，这个闭区间由$A$自身所决定。
如果$\vec{x}$碰巧是$A$的某个本征矢，则Reighleigh商即相应的本征值。
我们无法一步得到$A$的本征矢，那么退一步，我们寻求恰当的线性子空间，希望它们包含$\lambda_{min}(A), \lambda_{max}(A)$对应的本征矢。

假设现在有$k$个（$k\leq n$）正交归一的n维矢量$\{ \vec{q}_0, \vec{q}_1, \cdots \}$，定义
\begin{equation}
Q_k \equiv 
\begin{bmatrix}
\vec{q}_0, \vec{q}_1, \cdots, \vec{q}_{k-1}
\end{bmatrix},
\label{krylov-Qk}
\end{equation}
$Q_k$是一个$n \times k$的矩阵，因为$\vec{q}_i$的正交归一性，有$Q^\top_k Q_k = 1$。
$Q^\top_k A Q_k$是一个$k \times k$的方阵，那么，对任意非零$k \times 1$维矢量$\vec{y}$，有
\begin{equation}
r_{Q^\top_k A Q_k}(\vec{y}) 
= \frac{ \vec{y}^\top (Q^\top_k A Q_k) \vec{y} }{ \vec{y}^\top \vec{y} }
= \frac{ (Q_k \vec{y})^\top A (Q_k \vec{y}) }{ (Q_k \vec{y})^\top (Q_k \vec{y}) }
= r_A(Q_k \vec{y})
, ~~\vec{y} \neq 0
\end{equation}
注意，上式第二个等号使用了$\vec{Q}^\top_k Q_k = 1$
\footnote{下文中的Krylov子空间$\left\{ \vec{q}, A\vec{q}, \cdots, A^{k-1}\vec{q} \right\}$可以保证，$r_A(Q_k \vec{y})$确然越来越小，但在判断收敛条件时，在实际计算中我们需要对角化$Q^\top_k A Q_k$，要使$r_{Q^\top_k A Q_k}(\vec{y})$与$r_A(Q_k \vec{y})$，就需要$Q^\top Q = 1$。}。

根据上式与(\ref{eqn:Reighleigh-2}),
\begin{equation}
\lambda_{min} (Q^\top_k A Q_k) \leq r_{Q^\top_k A Q_k} (\vec{y}) = r_A(Q_k \vec{y}) <= \lambda_{max} (Q^\top_k A Q_k),
\end{equation}
$Q_k$的列矢$\vec{q}_0, \vec{q}_1, \cdots, \vec{q}_{k-1}$构成一个线性空间，是$n$维全空间中的一个子空间。
\begin{equation}
Q_k \vec{y} = \sum_i y_i \vec{q}_i,
\end{equation}
实则是$\vec{q}_i$的线性叠加，是$\vec{q}_i$构成的线性子空间中的一个点。
$\vec{x}=Q_k\vec{y}$时，$r_A(\vec{x})$在这个子空间里的极值，就是$\lambda_{min} (Q^\top_k A Q_k), \lambda_{max} (Q^\top_k A Q_k)$。

我们希望构造一系列子空间$Q_k, k=1,2,\cdots$，使得$r_A$在它们中的极值尽可能快地逼近$\lambda_{min}, \lambda_{max}$。
如果能做到这一点，我们就能很快地给出大型稀疏矩阵$A$的极端本征值的近似值。
最小的那些本征值可以是多体系统的低激发态，如果能得到它们的高质量近似，我们就可以摆宴席庆祝、感谢上天。

可以定义
\begin{eqnarray}
M_k &\equiv& \lambda_{max}( Q^\top_k A Q_k) = \max\limits_{\vec{y}\neq 0} r_{Q^\top_k A Q_k}(\vec{y}) = \max\limits_{\vec{y}\neq 0} r_A ( Q_k \vec{y}), \\
m_k &\equiv& \lambda_{min}( Q^\top_k A Q_k) = \min\limits_{\vec{y}\neq 0} r_{Q^\top_k A Q_k}(\vec{y}) = \min\limits_{\vec{y}\neq 0} r_A ( Q_k \vec{y}),
\end{eqnarray}

因为$Q_1, Q_2, \cdots, Q_{n}$的列矢量是逐渐增多的，它们的列矢量构成的线性空间是逐渐扩张的，所以有
\begin{eqnarray}
&& M_1 \leq M_2 \leq \cdots \leq M_n = \lambda_{max}(A), \\
&& m_1 \geq m_2 \geq \cdots \geq m_n = \lambda_{min}(A).
\end{eqnarray}

经过推导，可以得到
\begin{equation}
\nabla r_A(\vec{x}) = \frac{2}{\vec{x}^\top \vec{x}}( A \vec{x} - r_A(\vec{x})\vec{x} ).
\label{Rayleigh-gradient}
\end{equation}
所以，在Rayleigh商的极值点，有
\begin{equation}
\nabla r_A(\vec{x}) = 0,
\end{equation}
从而
\begin{equation}
A\vec{x} = r_A(\vec{x}) \vec{x}.
\end{equation}
所以，$r_A(\vec{x})$在$n$维线性空间中的极值点与A的本征值是一一对应的，本征值就是$r_A(\vec{x})$自己，这真是优美的性质！

如果选了一个初始的矢量$\vec{q}_1$，我们得到
\begin{eqnarray}
Q_1 &=& \begin{bmatrix}
\vec{q}_0
\end{bmatrix}, \\
M_1 &=& \max\limits_{y \neq 0} r_A ( Q_1 \vec{y}) = r_A (\vec{q}_0), \\
m_1 &=& \min\limits_{y \neq 0} r_A ( Q_1 \vec{y}) = r_A (\vec{q}_0).
\end{eqnarray}
我们想增加一个矢量$\vec{q}_2$，构造$Q_2$，使得$M_2 = \max\limits_{y \neq 0} r_A ( Q_2 \vec{y})$尽可能大，$m_2 = \min\limits_{y \neq 0} r_A ( Q_2 \vec{y})$尽可能小。
而$r_A(\vec{x})$的梯度为
\begin{equation}
\nabla r_A(\vec{x}) = \frac{2}{\vec{x}^\top \vec{x}}( A \vec{x} - r_A(\vec{x})\vec{x} )，
\end{equation}
所以我们可以设$\vec{q}_2 = A \vec{q}_1$，那么，$\{ \vec{q}_1, A\vec{q}_1 \}$张成的线性空间中，就包括了$\nabla r_A(\vec{q}_1)$，这意味着，这个线性空间中包含使得$r_A(\vec{x})$上升得最快的路径，也包括使它下降得最快的路径。
所以选择$\vec{q}_2 = A \vec{q}_1$是一个明智的决定，它可以使得$M_2$变得较大，$m_2$变得较小。

那么，如果$r_A(\vec{x})$在子空间$\{ \vec{q}_1, A\vec{q}_1 \}$中的极大值在$\vec{u} = u_1 \vec{q}_1 + u_2 A \vec{q}_1$处取得，那么下一个子空间最好包括
\begin{equation}
\nabla r_A(\vec{u}) = \frac{2}{\vec{u}^\top \vec{u}}( A \vec{u} - r(\vec{u})\vec{u} )
= \frac{2}{\vec{u}^\top \vec{u}}( u_1 A \vec{q}_1 + u_2 A^2 \vec{q}_1 - r(\vec{u})\vec{u} );
\end{equation}
如果$r_A(\vec{x})$在子空间$\{ \vec{q}_1, A\vec{q}_1 \}$中的极小值在$\vec{v} = v_1 \vec{q}_1 + v_2 A \vec{q}_1$处取得，那么下一个子空间最好包括
\begin{equation}
\nabla r_A(\vec{v}) = \frac{2}{\vec{v}^\top \vec{v}}( A \vec{v} - r(\vec{v})\vec{v} )
= \frac{2}{\vec{v}^\top \vec{v}}( v_1 A \vec{q}_1 + v_2 A^2 \vec{q}_1 - r(\vec{v})\vec{v} ).
\end{equation}
所以可以定
\begin{equation}
\vec{q}_3 = A^2 \vec{q}_1,
\end{equation}
就同时满足了上面两个要求，可以使得$M_3$尽量地大，$m_3$尽量地小。

类似地，定义Krylov子空间：
\begin{equation}
\kappa(A, \vec{q}_1, k) = \{ \vec{q}_1, A\vec{q}_1, \cdots, A^{k-1} \vec{q}_1 \},
\end{equation}
并用这些矢量按(\ref{krylov-Qk})构造$Q_k$.

这么做的话，$\kappa(A, \vec{q}_1, k)$中$r_A(\vec{x})$取最大值或者最小值处的梯度，都包括在$\kappa(A, \vec{q}_1, k)$子空间中。
所以使用Krylov子空间，可以使得$M_k$快速增大，$m_k$快速减小。

看起来，$\vec{q}_1$的选取，会对Krylov子空间的质量，即$M_k/m_k$增大/减小的速度，挺关键的。

\subsection{Lanczos三对角化}
上面定义了Krylov子空间，我们说明了，在给定$\vec{q}_0$的情况下，如果选取
\begin{equation}
\kappa(A, \vec{q}_0, k) = \{ \vec{q}_0, A\vec{q}_0, \cdots, A^{k-1} \vec{q}_0 \},
\end{equation}
可以让$Q^\top_k A Q_k$的极端本征值尽快地逼近$A$的极端本征值，其中$Q_k = [\vec{q}_0, \cdots, \vec{q}_k = A^{k-1}\vec{q}_0]$。
这是因为，增加$\vec{q}_k = A^{k-1}\vec{q}_0$以后，子空间$\kappa(A, \vec{q}_1, k)$就包括了
子空间$\kappa(A, \vec{q}_1, k-1)$中极值处的$\nabla r_A$方向。

这里面的一套操作与conjugate-gradient-method很相似，虽然我还没有理出关于相似性的清晰头绪来。

但是$\vec{q}_0, A\vec{q}_0, \cdots, A^{k-1} \vec{q}_0$是不正交不归一的，难以构造相似变换。
所以我们稍微改一下，设计一下$\vec{q}_k$，使得每次Krylov子空间扩张，既包括前一个子空间的极值处的$\nabla r_A$方向，同时也保证$\vec{q}_0, \cdots, \vec{q}_k, \cdots$正交归一。
这就是Lanczos方法，它将矩阵A或A的一部分三对角化。

定义$\vec{q}_{-1} = \beta_{-1} = 0$，给定归一化的$\vec{q}_0$，则$\forall k=0,\cdots,n-2$，
\begin{eqnarray}
\alpha_k &\equiv& \vec{q}^\top_k A \vec{q}_k, \\
\vec{r}_k &\equiv& (A-\alpha_kI)\vec{q}_k - \beta_{k-1} \vec{q}_{k-1}, \\
\beta_k &\equiv& \pm | \vec{r}_k |, \\
\vec{q}_{k+1} &\equiv& \vec{r}_k / \beta_k
\end{eqnarray}
因为$\vec{q}_{k+1}$包括了$A\vec{q}_k$的成分，所以$\kappa(A, \vec{q}_1, k)$就包括了
子空间$\kappa(A, \vec{q}_1, k-1)$中极值处的$\nabla r_A$方向。
$\vec{q}_{k+1}$也经过了归一化。

可以用归纳法证明，$\vec{q}_0, \vec{q}_1, \cdots$是正交的。
\begin{itemize}
\item[1] 首先，$\vec{r}_0 = A\vec{q}_0 - \alpha_0 \vec{q}_0$，所以有
\begin{equation}
\vec{q}^\top_0 \vec{r}_0 = \vec{q}^\top_0 A \vec{q}_0 - \alpha_0 = 0,
\end{equation}
即$\vec{q}_1, \vec{q}_0$正交。
而且，$\beta_{0} = \vec{q}^\top_1 \vec{r}_0 = \vec{q}^\top_1 A \vec{q}_0$。
\item[2] $\forall k = 1, \cdots$，如果$\beta_{k-1} = \vec{q}^\top_{k-1} \vec{q}_k$，而且
$\vec{q}_0, \cdots, \vec{q}_k$正交归一，那么用Lanczos方法做出来的$\vec{q}_{k+1}$和$\beta_{k}$满足以下几点
    \begin{itemize}
    \item[i] $\vec{q}_{k+1}$与$\vec{q}_k, \vec{q}_{k-1}$正交，因为
        \begin{equation}
        \vec{q}^\top_k \vec{r}_k = \alpha_k - \alpha_k = 0.
        \end{equation}
    \item[ii] $\beta_k = \vec{q}^\top_{k+1} A \vec{q}_k$，这是因为
        \begin{equation}
        \vec{q}^\top_{k+1} \vec{r}_k = \beta_k = \vec{q}^\top_{k+1} A \vec{q}_k.
        \end{equation}
    \item[iii] $\forall i=0, \cdots, k-2$，$\vec{q}_i$与$\vec{q}_{k+1}$正交，因为
        \begin{equation}
        \vec{q}^\top_i \vec{r}_k = \vec{q}^\top_i A \vec{q}_k = \vec{q}^\top_k A \vec{q}_i = \vec{q}^\top_k (\vec{r}_i + \alpha_i \vec{q}_i + \beta_{i-1}\vec{q}_{i-1})=0.
        \end{equation}
    \end{itemize}
\end{itemize}
所以，最终会有
\begin{equation}
AQ_n = Q_n T,
\end{equation}
其中$Q_n = [\vec{q}_1, \vec{q}_2, \cdots, \vec{q}_n]$，$T$为对称的三对角阵，对角元为$\alpha_k$，次对角元为$\beta_{k}$。
这就是Lanczos三对角化。

\paragraph{不变子空间}
如果在上面的三对角化过程中，发现
\begin{equation}
\vec{r}_k = (A - \alpha_k I)\vec{q}_k - \beta_{k-1}\vec{q}_{k-1} = 0,
\end{equation}
就说明遇到不变子空间了。
因为上式等价于
\begin{equation}
A \vec{q}_k = \alpha_k \vec{q}_k + \beta_{k-1} \vec{q}_{k-1},
\end{equation}
这意味着，A作用在$\kappa(A, \vec{q}_1, k)$上，得到的矢量还在这个子空间中，这就是不变子空间的定义了。

这样的话，三对角化就进行不下去。
当然，你可以把$Q_k$增补一些列矢量，扩充成$Q_n$，那么$Q_n^\top A Q_n$会成为分块矩阵，左上角那一块是三对角化的。

\subsection{Lanczos方法的实践}
我们把上面的$\vec{q}_0, \vec{q}_1, \cdots$叫做Lanczos矢量。
出于内存locality的考虑，我们不用在内存中储存所有的Lanczos矢量，而是将它们写到文件里。
在内存中，我们只需要储存两个矢量，让它们随着迭代往前滚动，并记录一路上的$\alpha_k, \beta_k$的值。
一边迭代，一边还可以对角化三对角矩阵$T_k$，如果所需极端本征值收敛了，就结束Lanczos过程，取出$T_k$的极端本征值。
如果需要相应的本征矢量，就读取文件中的Lanczos矢量，并且用$T_k$的本征矢作为权重，将$\vec{q}_0, \vec{q}_1, \cdots, \vec{q}_{k-1}$线性组合，就得到本征矢量。
如果极端本征值收敛了，就得到
\begin{equation}
\begin{bmatrix}
\vec{q}^\top_0 \\
\vdots \\
\vec{q}^\top_{k-1}
\end{bmatrix}
A
\begin{bmatrix}
\vec{q}_0 & \cdots & \vec{q}_{k-1}
\end{bmatrix}
= T_k
= V D V^\top,
\end{equation}
其中$V$的列矢量为$T_k$的本征矢，因为有
\begin{equation}
T_k V = V D,
\end{equation}
那么有
\begin{equation}
A
\begin{bmatrix}
\vec{q}_0 & \cdots & \vec{q}_{k-1}
\end{bmatrix}
V
= D \begin{bmatrix}
\vec{q}_0 & \cdots & \vec{q}_{k-1}
\end{bmatrix}
V,
\end{equation}
所以$\begin{bmatrix}
\vec{q}_0 & \cdots & \vec{q}_{k-1}
\end{bmatrix}
V$ 就构成A的近似本征矢。

具体算法可以总结如下：
\begin{shaded*}
\begin{itemize}
\item[1] 给定初始矢量$\vec{q}_0$，赋值给$\vec{\omega} = \vec{q}_0$，并储存$\vec{q}_0$，他们把它叫做“Lanczos pivot”。计算$\alpha_0 = \vec{q}^\top_0 A \vec{q}_0$。
\item[2] 计算$\vec{v} = A \vec{q}_0 - \alpha_0 \vec{q}_0$即前文中的$\vec{r}_0$， $\beta_0 = | \vec{v} | = |\vec{r}_0|$。
\item[3] 进入循环。第一次循环$k=1$。每次循环开始时，$\vec{v}=\vec{r}_{k-1}$, $\vec{\omega}=\vec{q}_{k-1}$。
    \begin{itemize}
    \item[i] 临时变量$\vec{t} = \vec{\omega}$；而$\vec{\omega} = \vec{v}/\beta_{k-1}$，即$\vec{q}_{k}$；$\vec{v} = - \beta_{k-1} \vec{t}$，即$-\beta_{k-1}\vec{q}_{k-1}$。
    \item[ii] 储存$\vec{q}_{k}$。
    \item[iii] $\vec{v} += A \vec{q}_{k}$。
    \item[iv] $\alpha_{k} = \vec{\omega}^\top \vec{v} = \vec{q}^\top_{k} A \vec{q}_{k}$
    \item[v] $\vec{v} -= \alpha_k \vec{q}_{k}$。到这里$\vec{v}$就是前面公式中的 $\vec{r}_{k}$。
    \item[vi] 正交化$\vec{v}$。因为数值误差，$\vec{v}$并不与前面的$\vec{q}_1, \cdots, \vec{q}_{k}$正交，所以需要做
    \begin{equation}
    \vec{v} -= \sum^{k}_{i=0} (\vec{v}^\top \vec{q}_i) \vec{q}_i,
    \end{equation}
    然后计算$\beta_{k} = | \vec{v} |$，即前面公式中的$|\vec{r}_{k}|$。
    \item[vii] $k=k+1$。
    \item[ix] $\alpha_0, \cdots, \alpha_{k-1}$作为对角元，$\beta_0, \cdots, \beta_{k-2}$作为次对角元构成的三对角阵，进行对角化，计算极端本征值与上次迭代的差别，如果差别小于精度要求，就结束循环，否则继续迭代。
    \end{itemize}
\end{itemize}
\end{shaded*}
根据Golub的书，Paige等人提出过更加有效的重新正交化方法，可以代替上面的3.vi，而不需要储存每一个$\vec{q}_i$在内存中。
在$n$比较大的时候，我这个比较简单直接的正交化方法因为memory locality的问题（因为储存每一个$\vec{q}_i$，可能会拖慢节奏。
下面是例程：
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#include<ctime>

#include"library.h"

double InnerProduct(int n, double *a, double *b){
        double y=0;
        for(int i=0;i<n;i++) y += a[i] * b[i];
        return y;
}

/*
 * Lanczos(...)         Lanczos method to diagonalize symmetric matrix A[n*n], the output is eigenval
ues in z[] and eigenvectors in vec[][]
 * int n                dimmension of the real symmetric matrix to be diagonalized
 * double *A            symmetric matrix
 * double *q0           initial pivot vector q0[n]
 * int nkeep            number of states to be kept
 * double *z            eigenvalues z[n] as output, whatever value it is when inputted.
 * double **vec         vec[nkeep][n] eigen vectors
 * double precision     precision
 */
void Lanczos(int n, double *A, int nkeep, double *z, double **vec, double precision){

        int i, j, k, l, max=nkeep * 100;
        double y, error = precision*10;
        double *omega = new double [n];
        double *v = new double [n];
        double *alpha = new double [n];
        double *beta = new double [n];
        double *d = new double [n];
        double *e = new double [n];

        for(i=0;i<n;i++) omega[i] = (double)rand()/RAND_MAX;// \vec{\omega} = \vec{q}_0
        y = sqrt( InnerProduct(n, omega, omega) );
        for(i=0;i<n;i++) omega[i] /= y;
        for(i=0;i<n;i++){
                y=0;
                for(j=0;j<n;j++){
                        y += A[i*n+j] * omega[j];// \vec{v} = A \vec{q}_0
                }
                v[i] = y;
        }
        alpha[0] = InnerProduct(n, omega, v);
        for(i=0;i<n;i++){
                v[i] -= alpha[0] * omega[i];// \vec{v} = \vec{r}_0 = A \vec{q}_0 - \alpha_0 \vec{q}_0
        }
        beta[0] = sqrt( InnerProduct(n, v, v) ); // beta[0] = | \vec{r}_0 |

        k = 1;

        double **q = new double *[max]; // Lanczos vectors q0, q1, ...
        for(i=0; i<max; i++) q[i] = new double [n];
        for(i=0; i<n; i++) q[0][i] = omega[i];// q[0] = q0
        // for test only
        /*
        double ** ove = new double *[n];
        for(i=0;i<n;i++){
                ove[i] = new double [n];
                for(j=0;j<n;j++) ove[i][j] = 0;
        } */
        double * eigket = new double [ max * max ];
        while(error > precision && k < max && k <= n ){
                cout<<"\n\n The "<<k<<"th Lanczos process:\n";
                for(i=0;i<n;i++){
                        y = omega[i];
                        omega[i] = v[i] / beta[k-1];// omega = \vec{q}_{k},
                        v[i] = - beta[k-1] * y;// v = - beta_{k-1} \vec{q}_{k-1}
                        q[k][i] = omega[i];// stores \vec{q}_k 
                }
                //-----------------------------------------------------------
                // for test only
                /*
                cout<<"ove:"<<endl;
                for(i=0;i<k+1;i++){
                        for(j=0;j<k+1;j++){
                                ove[i][j] = InnerProduct(n, q[i], q[j]);
                                if(fabs(ove[i][j])>1E-5) cout<<ove[i][j]<<",";
                                else cout<<"0,";
                        }
                        cout<<endl;
                }
                */
                //-----------------------------------------------------------

                for(i=0;i<n;i++){// O(n^2) is the most time consuming part
                        y=0;
                        for(j=0;j<n;j++){
                                y += A[i*n+j] * omega[j];
                        }
                        v[i] += y;// v = A q_{k+1} - beta_k q_k
                }
                y = 0;
                for(i=0;i<n;i++){
                        y += omega[i] * v[i];
                }
                alpha[k] = y;
                for(i=0;i<n;i++) v[i] -= alpha[k] * omega[i]; // v = A q_{k+1} - beta_k q_k - alpha_{k+1} q_{k+1} = r_{k+1}
                for(i=0;i<k+1;i++){// v - \sum_i (v^\top q_i) q_i // O(k*n)
                        y = InnerProduct(n, q[i], v);
                        for(j=0;j<n;j++)
                                v[j] -= y * q[i][j];
                }
                beta[k] = sqrt(InnerProduct(n,v,v)); // beta = | r_{k+1} |
                k = k+1;
                for(i=0;i<k;i++){
                        d[i] = alpha[i];
                        if(i==0) e[i] = 0;
                        else e[i] = beta[i-1];
                }
                /*
                cout<<"---------------------------------------------------------------"<<endl;
                cout<<"d:"; for(i=0;i<k;i++)cout<<d[i]<<"\t"; cout<<endl;
                cout<<"e:"; for(i=0;i<k;i++)cout<<e[i]<<"\t"; cout<<endl;
                */
                tqls(k, d, e, eigket);// O(k)
                sort_eig(k, eigket, d);
                //cout<<"tqls eig = "; for(i=0;i<k;i++)cout<<d[i]<<"\t"; cout<<endl;
                if( nkeep == n && k == n ){
                        cout<<"\t eigenvalues:"<<endl;
                        for(i=0;i<nkeep;i++) cout<<"\t "<<i+1<<"th eigenvalue: "<<d[i]<<endl;
                        break;
                }
                if( k > nkeep ){
                        y=0;
                        cout<<"\t eigenvalues:"<<endl;
                        for(i=0;i<nkeep;i++){
                                cout<<"\t "<<i+1<<"th eigenvalue: "<<d[i]<<endl;
                                y = (d[i] - z[i])*(d[i] - z[i]);
                        }
                        error = sqrt(y);
                        cout<<"\t error = "<<error<<" > "<<precision<<endl;
                }
                for(i=0;i<k;i++) z[i] = d[i];
        }
        if(k>max){
                cout<<"Lanczos: failed to converge after "<<max<<" iterations"<<endl;
                exit(1);
        }
        for(i=0;i<nkeep;i++){
                for(j=0;j<n;j++){
                        y=0;
                        for(l=0;l<k;l++){
                                y += eigket[i*k+l] * q[l][j];
                        }
                        vec[i][j] = y;
                }
        }
        delete [] beta;
        delete [] alpha;
        delete [] v;
        delete [] omega;
        delete [] d;
        delete [] e;
        delete [] eigket;
        for(i=0;i<max;i++) delete [] q[i];
        delete [] q;
}

extern "C" void dsyev_(char *a,char *b,int *n,double *A,int *nn,double *e,double *work,int *lwork,int *info);

int main(){

        int i,j,k,l, nkeep, n;
        double y;
        double *A, *a, *z, *q0, **vec;
        cout<<" matrix dimension = "; cin>>n;
        cout<<" nkeep = "; cin>>nkeep;

        srand(time(NULL));
        A = new double [n*n];
        for(i=0;i<n;i++){
                for(j=i;j<n;j++){
                        A[i*n+j] = (double)rand()/RAND_MAX;
                        A[j*n+i] = A[i*n+j];
                }
        }

        a = new double [n*n];
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        a[i*n+j] = A[i*n+j];
                }
        }
        z = new double [n]; for(i=0;i<n;i++)z[i]=0;

        q0 = new double [n];
        for(i=0;i<n;i++) q0[i] = (double)rand()/RAND_MAX;
        y = sqrt( InnerProduct(n, q0, q0) );
        for(i=0;i<n;i++) q0[i] /= y;

        vec = new double * [nkeep];
        for(i=0;i<nkeep;i++) vec[i] = new double [n];

        clock_t start = clock();
        Lanczos(n, A, nkeep, z, vec, 1E-6);
        clock_t end = clock();
        cout<<"Lanczos takes "<<(end-start)/CLOCKS_PER_SEC <<" s.\n";

        char alpha='v', beta='L';
        int lwork = 3 * n, info=0;
        double *e = new double [ n ];
        double *work = new double [ 3*n ];

        start = clock();
        dsyev_(&alpha, &beta, &n, a, &n, e, work, &lwork, &info);
        end = clock();
        cout<<"dsyev takes "<<(end-start)/CLOCKS_PER_SEC <<" s.\n";
        if(info<0){
                cout<<"The "<<abs(info)<<"th argument had an illegal value."<<endl;
                exit(1);
        }
        else if(info>0){
                cout<<"The algorithm failed to converge: "<<abs(info)<<" off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n";
                exit(1);
        }
        cout<<"lapack, eig = "; for(i=0;i<n;i++)cout<<e[i]<<"\t"; cout<<endl;

        // for test only
        /*
        for(i=0;i<nkeep;i++){
                cout<<"A v / lambda:";
                for(j=0;j<n;j++){
                        y=0;
                        for(k=0;k<n;k++){
                                y += A[j][k] * vec[i][k];
                        }
                        if( fabs(y/z[i]) > 1E-4 )cout<<y/z[i]<<",";
                }
                cout<<endl;
                cout<<"v:"; 
                for(j=0;j<n;j++){
                        if( fabs(vec[i][j]) > 1E-4 )cout<<vec[i][j]<<",";
                }       
                cout<<endl;
        }
        */

        for(i=0;i<nkeep;i++) delete [] vec[i]; delete [] vec;
        delete [] q0; delete [] z; delete [] a;
        delete [] A;

        return 0;
}
\end{lstlisting}
这个代码在linux下对比Lanczos方法与lapack包中的dsyev函数(三对角化＋对角化三对角阵)。
用到了library.h中的tqls与sort\_eig函数，在前文中有介绍。

上面的Lanczos函数中，还没有考虑invariant subspace的问题，即有可能$A^k \vec{q}_0 = 0$，那么会有$\vec{r}_k = 0$，所以需要新设一个垂直于$\vec{q}_0, \cdots, \vec{q}_{k-1}$的$\vec{q}_k$，重新启动Lanczos过程。

另外，因为计算中有误差，所以$\vec{q}_0, \cdots$的正交性很快会变差，所以Lanczos函数中，每次计算出$\vec{r}_{k+1}$，都重新正交化，即Ln102-106，
\begin{eqnarray}
\vec{r}_{k+1} -= \sum^k_{i=0}(\vec{r}^\top_{k+1} \vec{q}_i) \vec{q}_i.
\end{eqnarray}
Golub的书上有提到更好的重新正交化方法，应该是效率要更高。

\section{应用：求解1维无限深势井}
在第八章线性方程组的开头，我们提到过，利用有限差分，二阶常微分方程可以转化为线性方程组求解问题。
那么，如果微分方程中有待定的特征值，则可以考虑将微分方程转化为矩阵本征值问题。
1维无限深势井，就是一个这样的例子。
\begin{equation}
V(x) = \left\{
\begin{aligned}
&0, ~~~ 0<x<a\\
&\infty, ~~else
\end{aligned}
\right.
\end{equation}
在$[0,a]$区间内，薛定谔方程为
\begin{equation}
- \frac{ \hbar^2 }{2m} \frac{ d^2 }{dx^2} \psi = E \psi，
\end{equation}
相应的本征值和解为
\begin{eqnarray}
E_n &=& \frac{ n^2 \pi^2 \hbar^2 }{ 2 m a^2}, \\
\psi_n &=& \sqrt{\frac{2}{a}} \sin (n \pi x /a  ).
\end{eqnarray}

在数值中，为了简便，我们取$\frac{\hbar^2}{2m} = a = 1$，则薛定谔方程变为
\begin{equation}
- \frac{d^2}{dx^2} \psi = E \psi, \label{linear-1d-infinite-well}
\end{equation}
其解为
\begin{eqnarray}
E_n &=& n^2 \pi^2 , \\
\psi_n &=& \sqrt{2} \sin (n \pi x  ).
\end{eqnarray}
我们可以将$[0,a]$，即$[0,1]$区间等分$n+1$份，则有$n$个节点$x_1, \cdots, x_{n}$，当然，还有$x_{0} = 0, x_{n+1} = 1$。
相应地，记各个节点上$\psi$的数值近似为$y_{0}, y_1, ..., y_{n+1}$，由边界条件，我们有
\begin{equation}
y_{0} = 0, ~~ y_{n+1} = 0.
\end{equation}
利用有限差分公式
\begin{eqnarray}
\frac{d^2}{dx^2} \psi (x) |_{x_i} &\approx& ( \psi (x_{i+1}) + \psi(x_{i-1}) - 2\psi(x_i) )/h^2, 
\end{eqnarray}
得到$y_0, \cdots, y_{n+1}$之间的关系。
在$x_1, x_2, \cdots, x_n$各个节点上，(\ref{linear-1d-infinite-well})转化为
\begin{eqnarray}
-\frac{1}{h^2}y_0 + \frac{2}{h^2}y_1 - \frac{1}{h^2}y_2 &=& E y_1, \\
-\frac{1}{h^2}y_1 + \frac{2}{h^2}y_2 - \frac{1}{h^2}y_3 &=& E y_2, \\
-\frac{1}{h^2}y_2 + \frac{2}{h^2}y_3 - \frac{1}{h^2}y_4 &=& E y_3, \\
\cdots&& \\
-\frac{1}{h^2}y_{n-1} + \frac{2}{h^2}y_n - \frac{1}{h^2}y_{n+1} &=& E y_n, \\
\end{eqnarray}
将$y_0 = y_{n+1} = 0$代入其中，就得到一个矩阵形式的本征值问题：
\begin{equation}
\begin{bmatrix}
\frac{2}{h^2}  & - \frac{1}{h^2}  \\
-\frac{1}{h^2} & \frac{2}{h^2}   & - \frac{1}{h^2}\\
				& -\frac{1}{h^2} & \frac{2}{h^2}   & - \frac{1}{h^2} \\
				&				 & \vdots          & \vdots        & \vdots \\
				&				&				&	- \frac{1}{h^2} &  	\frac{2}{h^2}	
\end{bmatrix}
\begin{bmatrix}
y_1 \\
y_2 \\
y_3 \\
\vdots \\
y_n
\end{bmatrix}
= E
\begin{bmatrix}
y_1 \\
y_2 \\
y_3 \\
\vdots \\
y_n
\end{bmatrix}
\end{equation}
对角化上式左侧的矩阵，就得到了能量和波函数。

下面的函数，即实现上面的算法，
\begin{lstlisting}
/* ISW1d:       解出1维无限深方势井本征值和波函数。假定 hbar^2/2m = a = 1，ISW1d: Infinite-Square-Well-1-dimension。解出的波函数写进"ISW1d_wave.txt"中。
 * int n:       (0,1)上的节点数
 * double *e:   本征值，从小到大排列
 */
double ISW1d(int n, double *e){//1d Infinite Square Well

        int i,j;
        double h = 1.0/(n+1);// 将区域 [0,1] 分为 n+1 份，所以其中有 n 个节点, 步长为h
        double *a = new double [n*n];//课件中的矩阵
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        if(i==j){
                                a[i*n+j] = 2/h/h;
                        }
                        else if( abs(i-j)==1 ) {
                                a[i*n+j] = -1/h/h;
                        }
                        else a[i*n+j]=0;
                }
        }
        
        jacobi(n, a, e);//对角化
        sort_eig(n, a, e);//按本征值从小到大，将本征值和本征矢量排序

        // 将波函数输出到ISW1d_wave.txt
        ofstream fp("ISW1d_wave.txt");
        for(i=0;i<n+1;i++)fp<<"0.0\t"; fp<<endl;
        for(i=0;i<n;i++){
                fp<<setprecision(4)<<(i+1)*h<<"\t";
                for(j=0;j<n;j++)
                        fp<<setprecision(4)<<a[j*n+i]<<"\t";
                fp<<endl;
        }
        fp<<"1.0\t";
        for(i=0;i<n;i++)fp<<"0.0\t"; fp<<endl;
        fp.close();

        delete [] a;//释放动态内存
        
        return 0;
}
\end{lstlisting}
这个函数可以求出本征值来。
如果要得到更精确的结果，我们需要增大$n$，减小步长$h$。
但是计算量随着$n$的变大，迅速地变大，是否有更好的办法？

我们可以这样看待这个问题。
求得的$E_i$都是步长$h$的函数，当$h$趋于0的时候，$E_i$一定会趋于理论值，只是我们算不动太大的矩阵。
那么，我们可以选取不同的步长值，$h, h/2, h/4$，分别得到$E^h_i, E^{h/2}_i, E^{h/4}_i$，然后做多项式外推，推测$h=0$时$E_i$的值，也会提高我们的精度。

在下面的代码中，我们调用第二部分写的Nevill函数，用Nevill外推法推得$E_i$的值。
\begin{lstlisting}
int main(){
        
        int n=30, i, j;
        double y, dy;

        double x_data[3], eig[10][3];//用作Nevill外推的输入量

        for(i=0;i<3;i++){
                x_data[i] = 1.0/(n+1);//步长
                double *e = new double [n];
                ISW1d(n, e);
                cout<<"e:"; for(j=0;j<10;j++)cout<<e[j]/M_PI/M_PI<<","; cout<<endl;
                for(j=0;j<10;j++) eig[j][i] = e[j];//本征值
                delete [] e;
                n*=2;
        }

        cout<<"Extrapolation:"<<endl;
        for(i=0;i<10;i++){
                Nevill(3, x_data, eig[i], 0.0, y, dy);// y为外推得到的 E_{i+1}
                cout<<"E"<<(i+1)<<"/pi/pi = "<<y/M_PI/M_PI<<endl;
        }
        return 0;
}
\end{lstlisting}
外推得到的前10个能级如下
\begin{table}[h]
\centering
\begin{tabular}{ccc}
\hline\hline
i & $E_i/\pi^2$ & 理论值 \\
\hline
1 & 1.000		& 1 \\
2 & 4.000		& 4 \\
3 & 9.000		& 9 \\
4 & 16.000		& 16 \\
5 & 25.001		& 25 \\
6 & 36.003		& 36 \\
7 & 49.008		& 49 \\
8 & 64.017		& 64 \\
9 & 81.035		& 81 \\
10 & 100.066 	& 100 \\
\hline\hline
\end{tabular}
\caption{一维无限深势阱前10个能级，保留三位小数。利用有限差分，将薛定谔方恒转化为矩阵本征值问题，然后用不同的步长计算$E_i$，再外推得到步长$h\rightarrow 0$时的$E_i$值。}
\end{table}
上面的计算中，用到的矩阵的最大维度为$120 \times 120$，但已经推得非常精确的结果。
下图中是$n=120$时求得的最低的几个能级的波函数
\begin{figure*}[h]
\centering
\includegraphics[width=0.8\textwidth]{ISW1d}
\caption{一维无限深势井的波函数。}
\end{figure*}

\section{第3部分作业}
1. 怎样用高斯消元法求矩阵的逆矩阵？是否可以求本征值？（简答，将回答写在邮件正文即可）\\
2. 利用ludcmp函数，求矩阵的行列式：
\begin{equation}
A=
\begin{pmatrix}
1   & 1 &   1   &   1\\
1   & 2 &   2   &   2\\
1   & 2 &   3   &   3\\
1   & 2 &   3   &   4
\end{pmatrix}.
\end{equation}
*3. 利用ludcmp函数，求矩阵的逆矩阵：
\begin{equation}
A=
\begin{pmatrix}
1   & 1 &   1   &   1\\
1   & 2 &   2   &   2\\
1   & 2 &   3   &   3\\
1   & 2 &   3   &   4
\end{pmatrix}.
\end{equation}
4. 实验室制备了100个完全相同的放射性X原子，它们在单位时间内衰变的概率相同。每隔十分钟数一遍剩下的 X 原子数,得到数据：
(9:00, 100), (9:10, 71), (9:20, 50), (9:30, 37), 
(9:40, 26), (9:50, 17), (10:00, 11), (10:10, 9)
试用最小二乘法拟合这种原子的半衰期。

**5. 试编写jacobi迭代法程序，求解以下方程组的解
\begin{equation}
\begin{pmatrix}
1   & 1 &   1   &   1\\
1   & 2 &   2   &   2\\
1   & 2 &   3   &   3\\
1   & 2 &   3   &   4
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
x_4
\end{pmatrix}
=\begin{pmatrix}
4\\
7\\
9\\
10
\end{pmatrix}.
\end{equation}

6. 利用jacobi方法，求以下矩阵的本征值与本征矢量：
\begin{equation}
A = \begin{pmatrix}
0.575155 & 0.878075 & 0.939033\\
0.878075 & 0.445565 & 0.99726\\
0.939033 & 0.99726 & 0.957276
\end{pmatrix}
\end{equation}
% 答案：eigenvalues: -0.397167,-0.183899,2.55906,
%eigenkets: 
%-0.482712,0.830803,-0.277047,
%0.687128,0.163144,-0.707983,
%0.542996,0.532118,0.649619,

7. 试用课件中展示的有限差分+矩阵对角化+Nevill外推方法，求解1维谐振子势
\begin{equation}
V(x) = \frac{1}{2} k x^2,
\end{equation}
为了叙述简便，不妨取$\hbar = m = k = 1$, 则薛定谔方程为
\begin{eqnarray}
- \frac{d^2}{dx^2} \psi + x^2 \psi &=& 2 E \psi, \\
                                 E &=& (n+\frac{1}{2}).
\end{eqnarray}
我们可以选取区间$[-5,5]$，将区间等分为$n+1$份，并假设$\psi(-5) = \psi(5) = 0$，然后进行类似于讲义中类似于一维无限深势井的处理，\\
1)试推得前3个本征值，\\
*2)画出前三个本征态的波函数，并做图。

